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The double pulsar system, PSR J0737-3039A/B, is unique in that both neutron 
stars are detectable as radio pulsars. This, combined with significantly higher 
mean orbital velocities and accelerations when compared to other binary pul- 
sars, suggested that the system would become the best available testbed for 
general relativity and alternative theories of gravity in the strong-field regime. 
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Here we report on precision timing observations taken over the 2.5 years since 
its discovery and present four independent strong-field tests of general rela- 
tivity. Use of the theory-independent mass ratio of the two stars makes these 
tests uniquely different from earlier studies. By measuring relativistic correc- 
tions to the Keplerian discription of the orbital motion, we find that the "post- 
Keplerian" parameter s agrees with the value predicted by Einstein's theory 
of general relativity within an uncertainty of 0.05%, the most precise test yet 
obtained. We also show that the transverse velocity of the system's center of 
mass is extremely small. Combined with the system's location near the Sun, 
this result suggests that future tests of gravitational theories with the double 
pulsar will supersede the best current Solar-system tests. It also implies that 
the second-born pulsar may have formed differently to the usually assumed 
core-collapse of a helium star. 

Introduction. Einstein's general theory of relativity (GR) has so far passed all experimental 
tests with flying colours O, with the most precise tests achieved in the weak-field gravity 
conditions of the Solar System (0131. However, it is conceivable that GR breaks down under 
extreme conditions such as strong gravitational fields where other theories of gravity may apply 
0. Predictions of gravitational radiation and self-gravitational effects can only be tested using 
massive and compact astronomical objects such as neutron stars and black holes. Studies of 
the double-neutron- star binary systems, PSR B1913-I-16 and PSR B 1534-1- 12, have provided the 
best such tests so far, confirming GR at the 0.2% and 0.7% level, respectively (|5l|6|) ^ The 

recently discovered double pulsar system, PSR J0737-3039A/B, has significantly higher mean 

'Stairs et al. (2002, ref. (|6j) find an agreement of their measured values for PSR B1534+12 with GR at the 
0.05% level, but the measurement uncertainty on the most precisely measured parameter in the test, s, is only 
0.7%. 
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orbital velocities and accelerations than either PSR B1913+16 or PSR B1534+12 and is unique 
in that both neutron stars are detectable as radio pulsars 

PSR J0737-3037A/B consists of a 22-ms period pulsar, PSR J0737-3039A (henceforth 
called A), in a 2.4-hr orbit with a younger 2.7-s period pulsar, PSR J0737— 3039B (B). Soon 
after the discovery of A Q, it was recognised that the orbit's orientation, measured as the 
longitude of periastron cu, was changing in tine with a very large rate of cj = dcu/dt ~ 17° 
yr \ which is four times the corresponding value for the Hulse-Taylor binary, PSR B 19 13+16 
This immediately suggested that the system consists of two neutron stars, a conclusion 
confirmed by the discovery of pulsations from B ©. The pulsed radio emission from B has 
a strong orbital modulation, both in intensity and in pulse shape. It appears as a strong radio 
source only for two intervals, each of about 10-min duration, while its pulsed emission is rather 
weak or even undetectable for most of the remainder of the orbit (0111). 

In double-neutron-star systems, especially those having short orbital periods, observed pulse 
arrival times are significantly modified by relativistic effects which can be modelled in a theory- 
independent way using the so-called "Post-Keplerian" (PK) parameters (7^. These PK param- 
eters are phenomenological corrections and additions to the simple Keplerian description of the 
binary motion, describing for instance a temporal change in period or orientation of the orbit, or 
an additional "Shapiro-delay" that occurs due to the curvature of space-time when pulses pass 
near the massive companion. The PK parameters take different forms in different theories of 
gravity and so their measurement can be used to test these theories dTTllTl) . For point masses with 
negligible spin contributions, GR predicts values for the PK parameters which depend only on 
the two a priori unknown neutron-star masses and the precisely measurable Keplerian parame- 
ters. Therefore measurement of three (or more) PK parameters provides one (or more) tests of 
the predictive power of GR. For the double pulsar we can also measure the mass ratio of the 
two stars, R = nip^/m-Q = xb/xa- The ability to measure this quantity provides an important 
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constraint because in GR and other theories this simple relationship between the masses and 
semi-major axes is valid to at least first post-Newtonian (IPN) or (f /c)^ order (\12\\liy 

Observations. Timing observations of PSR J0737— 3039A/B have been undertaken using the 
64-m Parkes radio telescope in New South Wales, Australia, the 76-m Lovell radio telescope 
at Jodrell Bank Observatory (JBO), UK, and the 100-m Green Bank Telescope (GBT) in West 
Virginia, USA, between 2003 April and 2006 January. 

At Parkes, observations were carried out in bands centred at 680 MHz, 1374 MHz and 
3030 MHz. While timing observations were frequent after the discovery of the system, later 
observations at Parkes were typically conducted every 3-4 weeks, usually covering two full 
orbits per session. Observations at the GBT were conducted at monthly intervals, with each 
session consisting of a 5- to 8-hour track (i.e., 2 to 3 orbits of the double pulsar). Typically, the 
observing frequencies were 820 and 1400 MHz for alternate sessions. Occasionally, we also 
performed observations at 340 MHz, in conjunction with pulse profile studies to be reported 
elsewhere. In addition, we conducted concentrated campaigns of five 8-hour observing sessions, 
all at 820 MHz, in 2005 May and 2005 November. Observations at JBO employed the 76-m 
Lovell telescope. Most data were recorded at 1396 MHz, while some observing sessions were 
carried out at the lower frequency of 610 MHz. The timing data obtained at Jodrell Bank 
represent the most densely sampled dataset but, because of the limited bandwidth, requiring 
longer integration times per timing point. The Parkes dataset is the longest one available and 
hence provides an excellent basis for investigation of secular timing terms. 

The time-series data of all systems were folded modulo the predicted topocentric pulse 
period. The adopted integration times were 30 s for pulsar A (180 s for JBO data) and 300 s for 
pulsar B. For A, these integration times reflect a compromise between producing pulse profiles 
with adequate signal-to-noise ratio and sufficient sampling of the orbit to detect and resolve 
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phenomena that depend on orbital phase, such as the Shapiro delay. The integration time for B 
corresponds to about 108 pulse periods and is a compromise between the need to form a stable 
pulse profile while resolving the systematic changes seen as a function of orbital phase. 

Timing measurements. For each of the final profiles, pulse times-of-arrival (TOAs) were 
computed by correlating the observed pulse profiles with synthetic noise-free templates (see 
Fig. 1 in (I7?t . cf. ref. Q). A total of 131,416 pulse TOAs were measured for A while 507 
TOAs were obtained for B. For A, the same template was used for all observations in a given 
frequency band, but different templates were used for widely separated bands. We note that our 
observations still provide no good evidence for secular evolution of A's profile (175b despite the 
predictions of geodetic precession. The best timing precision was obtained at 820 MHz with 
GASP backend (see ref. Hlh for details of this and other observing systems) on the GBT, with 
typical TOA measurement uncertainties for pulsar A of 18 /is for a 30-s integration. 

For B, because of the orbital and secular dependence of its pulse profile (9), different tem- 
plates were also used for different orbital phases and different epochs. A matrix of B templates 
was constructed, dividing the data set into 3-month intervals in epoch and 5-minute intervals in 
orbital phase. The results for the 29 orbital phase bins were studied, and it was noticed that, 
while the profile changes dramatically and quickly during the two prominent bright phases, the 
profile shape is simpler and more stable at orbital phases when the pulsar is weak. This appar- 
ent stability at some orbital phases cannot be attributed to a low signal-to-noise ratio as secular 
variations in the pulse shape are still evident. Consequently, the orbital phase was divided into 
five groups of different lengths to which the same template (for a given 3-month interval) was 
applied as shown in Fig. 2 of dT?!) . In the final timing analysis, data from the two groups repre- 
senting the bright phases (IV & V in Fig. 2 of dTib ) were excluded to minimize the systematic 
errors caused by the orbital profile changes. Also, because of signal-to-noise and radio inter- 
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ference considerations, only data from Parkes and the GBT BCPM backend were used in the B 
timing analysis. 

All TOAs were transferred to Universal Coordinated Time (UTC) using the Global Posi- 
tional System (GPS) to measure offsets of station clocks from national standards and Circu- 
lar T of the BIPM to give offsets from UTC, and then to the nominally uniform Terrestrial 
Time TT(BIPM) timescale. These final TOAs were analysed using the standard software pack- 
age TEMPO (I73t . fitting parameters according to the relativistic and theory independent timing 
model of Damour & Deruelle ilTWlOb . In addition to the DD model, we also applied the "DD- 
Shapiro" (DDS) model introduced by Kramer et al. (ref. (7^)). The DDS model is a modification 
of the DD model designed for highly inclined orbits. Rather than fitting for the Shapiro param- 
eter s, the model uses the parameter Zg = — ln(l — s) which gives a more reliable determination 
of the uncertainties in Zg and hence in s. We quote the final result for the more commonly used 
parameter s and note that its value computed from Zg is in good agreement with the value ob- 
tained from a direct fit for s within the DD model. Derived pulsar and binary system parameters 
are listed in Table 1 . 

In the timing analysis for pulsar B, we used an unweighted fit to avoid biasing the fit toward 
bright orbital phases. Uncertainties in the timing parameters were estimated using Monte Carlo 
simulations of fake data sets for a range of TOA uncertainties, ranging from the minimum esti- 
mated TOA error to its maximum observed value of about 4 ms. For B, we also fitted for offsets 
between datasets derived from different templates in the fit since the observed profile changes 
prevent the establishment of a reliable phase relationship between the derived templates. This 
precludes a coherent fit across the whole orbit and hence limits the final timing precision for B. 
It cannot yet be excluded that different parts of B's magnetosphere are active and responsible 
for the observed emission at different orbital phases. 

In the final fit, we adopted the astrometric parameters and the dispersion measure derived 
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for A and held these fixed during the fit, since A's shorter period and more stable profile give 
much better timing precision than is achievable for B. Except for the semi-major axis which 
is only observable as the projection onto the plane-of-the-sky xb — (ob/c) sini, where i is the 
orbital inclination angle, we also adopted A's Keplerian parameters (with 180° added to loa) and 
kept these fixed. We also adopted the PK parameter to from the A fit since logically this must 
be identical for the two pulsars; this equality therefore does not implicitly make assumptions 
about the validity of any particular theory of gravity (see next section). The same applies for 
Pb. In contrast, the PK parameters 7, s and r are asymmetric in the masses and their values 
and interpretations differ for A and B. In practical terms, the relatively low timing precision 
for B does not require the inclusion of 7, s, r or Pb in the timing model. We can however 
independently measure cub, obtaining a value of 16.96 ± 0.05 deg yr"^, consistent with the 
more accurately determined value for A. 

Since the overall precision of our tests of GR is currently limited by our ability to measure 
Xb and hence the mass ratio R = mx/mB = xb/x^ (see below), we adopted the following 
strategy to obtain the best possible accuracy for this parameter. We used the whole TOA data 
set for B in order to measure B's spin parameters P and P, given in Table 1. These parameters 
were then kept fixed for a separate analysis of the concentrated 5 -day GBT observing sessions 
at 820 MHz. On the timescale of the long-term profile evolution of B, each 5-day session 
represents a single-epoch experiment and hence requires only a single set of profile templates. 
The value of xb obtained from a fit of this parameter only to the two 5-day sessions is presented 
in Table 1. 

Because of the possible presence of unmodelled intrinsic pulsar timing noise and because 
not all TOA uncertainties are well understood, we adopt the common and conservative pulsar- 
timing practice of reporting twice the parameter uncertainties given by TEMPO as estimates of 
the 1-(T uncertainties. While we believe that our real measurement uncertainties are actually 
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somewhat smaller than quoted, this practice facilitates the comparison with previous tests of 
GR using pulsars. The timing model also includes timing offsets between the datasets for the 
different instruments represented by the entries in Table 1 in {13). The final weighted rms 
post-fit residual is 54.2/iS. In addition to the spin and astrometric parameters, the Keplerian 
parameters of A's orbit and five PK parameters, we also quote a tentative detection of a timing 
annual parallax which is consistent with the dispersion-derived distance. Further details are 
given in ref. (I7?t . 

Tests of general relativity. Previous observations of PSR J0737— 3039A/B (001 resulted in 
the measurement of R and four PK parameters: the rate of periastron advance u, the grav- 
itational redshift and time dilation parameter 7, and the Shapiro-delay parameters r and s. 
Compared to these earlier results, the measurement precision for these parameters from PSR 
J0737— 3039A/B has increased by up to two orders of magnitude. Also, we have now mea- 
sured the orbital decay, P^. Its value, measured at the 1.4% level after only 2.5 years of timing, 
corresponds to a shrinkage of the pulsars' separation at a rate of 7mm per day. Therefore, we 
have measured five PK parameters for the system in total. Together with the mass ratio R, we 
have six different relationships that connect the two unknown masses for A and B with the ob- 
servations. Solving for the two masses using R and a one PK parameter, we can then use each 
further PK parameter to compare its observed value with that predicted by GR for the given 
two masses, providing four independent tests of GR. Equivalently, one can display these tests 
elegantly in a "mass-mass" diagram (Fig. 1). Measurement of the PK parameters gives curves 
on this diagram that are in general different for different theories of gravity but which should 
intersect in a single point, i.e., at a pair of mass values, if the theory is valid (llL 

As shown in Fig. 1, we find that all measured constraints are consistent with GR. The 
most precisely measured PK parameter currently available is the precession of the longitude 
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of periastron, uj. We can combine this with the theory-independent mass ratio R to derive 
the masses given by the intersection region of their curves: rriA = 1.3381 ± 0.0007 M0 
and me = 1.2489 ± 0.0007 Mq.^ Table 2 lists the resulting four independent tests that are 
currently available. All of them rely on comparison of our measured values of s, r, 7 and 
Pb with predicted values based on the masses defined by the intersection of the allowed re- 
gions for iu and R in the niA-rnB plane. The calculation of the predicted values is some- 
what complicated by the fact that the orbit is nearly edge-on to the line of sight, so that the 
formal intersection region actually includes parts of the plane disallowed by the Keplerian 
mass functions of both pulsars (see Fig. 1). To derive legitimate predictions for the various 
parameters, we used the following Monte Carlo method. A pair of trial values for iu and 
xb (and hence R and the B mass function) is selected from gaussian distributions based on 
the measured central values and uncertainties. (The uncertainty on xa is very small and is 
neglected in this procedure.) This pair of trial values is used to derive trial masses niA and 
thb, using the GR equation Co = 3{^)~^/^{TqM)'^/^ (1 - e^)"\ where M = tua + and 
Tq = GMq/c^ = 4.925490947 /iS, and the mass-ratio equation mx/m-Q = xb/xa- If this trial 
mass pair falls in either of the two disallowed regions (based on the trial mass function for B) 
it is discarded. This procedure allows for the substantial uncertainty in the B mass function. 
Allowed mass pairs are then used to compute the other PK parameters, assuming GR. This pro- 
cedure is repeated until large numbers of successful trials have accumulated. Histograms of the 
PK predictions are used to compute the expectation value and 68% confidence ranges for each 
of the parameters. These are the values given in Table 2. 

The Shapiro delay shape illustrated in Fig. 2 gives the most precise test, with Sobs/spred = 

^The true masses will deviate from these values by an unknown, but essentially constant, Doppler factor, 
probably of order 10^"^ or less (iOi. Moreover, what is measured is a product containing Newton's gravitational 
constant G. The relative uncertainty of G of 1.5 x 10^"' limits our knowledge of any astronomical mass in kilograms 
but since the product Tq = GMq jc" = 4.925490947/is is known to very high precision, masses can be measured 
precisely in solar units. 
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0.99987 ± 0.00050.^ This is by far the best available test of GR in the strong-field limit, having 

a higher precision than the test based on the observed orbit decay in the PSR B1913+16 system 

with a 30-year data span As for the PSR B1534+12 system the PSR J0737-3039A/B 

Shapiro-delay test is complementary to that of B 1913-1-16 since it is not based on predictions 

relating to emission of gravitational radiation from the system (l2Qt . Most importantly, the four 

tests of GR presented here are qualitatively different from all previous tests because they include 

one constraint (R) that is independent of the assumed theory of gravity at the IPN order. As a 

result, for any theory of gravity, the intersection point is expected to lie on the mass ratio line in 

Fig. 1 . GR also passes this additional constraint. 

Li estimating the final uncertainty of xb and hence of R, we have considered that geodetic 

precession will lead to changes to the system geometry and hence changes to the aberration of 

the rotating pulsar beam. The effects of aberration on pulsar timing are usually not separately 

measurable but are absorbed into a redefinition of the Keplerian parameters. As a result, the 

observed projected sizes of the semi-major axes, x^^, differ from the intrinsic sizes, a;™^ by 

a factor (1 + g). The quantity ca depends for each pulsar A and B on the orbital period, 

the spin frequency, the orientation of the pulsar spin and the system geometry (|77Tl . While 

aberration should eventually become detectable in the timing, allowing the determination of a 

further PK parameter, at present it leads to an undetermined deviation of x"^^ from x™*, where 

the latter is the relevant quantity for the mass ratio. The parameter g scales with pulse period 

and is therefore expected to be two orders of magnitude smaller for A than for B. However, 

because of the high precision of the A timing parameters, the derived value x'p^^ may already 

be significantly affected by aberration. This has (as yet) no consequences for the mass ratio 

R = x^^/x'^^, as the uncertainty in R is dominated by the much less precise x^^. We can 

explore the likely aberration corrections to x"^^ for various possible geometries. Using a range 
^Note, s has the same relative uncertainty as our determination of the masses. 
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of values given by studies of the double pulsar's emission properties (U7T) . we estimate ~ 
10"^ and eg ~ 10^"^. The contribution of aberration therefore is at least one order of magnitude 
smaller than our current timing precision. In the future this effect may become important, 
possibly limiting the usefulness of R for tests of GR. If the geometry cannot be independently 
determined, we could use the observed deviations of R from the value expected within GR to 
determine and hence the geometry of B. 

Space motion and inclination of the orbit. Because the measured uncertainty in de- 
creases approximately as T^^^, where T is the data span, we expect to improve our test of 
the radiative aspect of the system to the 0.1% level or better in about five years' time. For the 
PSR B1913+16 and PSR B1534+12 systems, the precision of the GR test based on the orbit- 
decay rate is severely limited both by the uncertainty in the differential acceleration of the Sun 
and the binary system in the Galactic gravitational potential and the uncertainty in pulsar dis- 
tance ([221 16)). For PSR J0737— 3039A/B, both of these corrections are very much smaller than 
for these other systems. Based on the measured dispersion measure and a model for the Galactic 
electron distribution PSR J0737— 3039A/B is estimated to be about 500 pc from the Earth. 
From the timing data we have measured a marginally significant value for the annual parallax, 
3 ± 2 mas, corresponding to a distance of 200 — 1000 pc (Table 1), which is consistent with the 
dispersion-based distance that was also used for studies of detection rates in gravitational wave 
detectors Q. The observed proper motion of the system (Table 1) and differential acceleration 
in the Galactic potential then imply a kinematic correction to Pb at the 0.02% level or less. 
Independent distance estimates also can be expected from measurements of the annual parallax 
by Very Long Baseline Interferometry (VLBI) observations, allowing a secure compensation 
for this already small effect. A measurement of Pb at the 0.02% level or better will provide 
stringent tests for alternative theories of gravity. For example, limits on some scalar-tensor 
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theories will surpass the best current Solar-system tests (05t . 

In GR, the parameter s can be identified with sin i where i is the inclination angle of the 
orbit. The value of s given in Table 1 corresponds to i = 88°.69loo fg. Based on scintillation 
observations of both pulsars over the short time interval when A is close to superior conjunction, 
Coles et al. S26[ derived a value for \i — 90° | of 0°.29 ± 0°.14. This is consistent with our 
measurement only at the 3-cr level. As mentioned above, we used the DDS model to solve for 
the Shapiro delay. Fig. 3 shows the resulting contours in the Zg - tub plane. The value and 
uncertainty range for s quoted in Table 1 correspond to the peak and range of the 68% contour. 
Because of the non-linear relationship between Zg and s, the uncertainty distribution in s (and 
hence in i) corresponding to these contours is very asymmetric with a very steep edge on the 90° 
side. Only close to the 99% confidence limit is the timing result consistent with the scintillation- 
derived value of \i — 90° | of 0°.29 ± 0°.14 i,26). We note that the scintillation measurement is 
based on the correlation of the scintillation fluctuations of A and B over the short interval when 
A is close to superior conjunction (i.e., behind B). In contrast, the measurement of i from timing 
measurements depends on the detection of significant structure in the post-fit residuals after a 
portion of the Shapiro delay is absorbed in the fit for xa (E3- As shown in Fig. 2, the Shapiro 
delay has a signature that is spread over the whole orbit and hence can be cleanly isolated. We 
also examined the effects on the Shapiro delay of using only low- or high-frequency data, and 
found values of s consistent withing the errors in each case. The scintillation result is based on 
the plasma properties of the interstellar medium and may also be affected by possible refraction 
effects in B's magneto sphere. We believe that the timing result is much less susceptible to 
systematic errors and is therefore more secure. 

Scintillation observations have also been used to deduce the system transverse velocity. 
Ransom et al. (l2Sb derive a value of 141 ± 8.5 km s^^ while Coles et al. {26) obtain 66 ± 15 km 
s^^ after considering the effect of anisotropy in the scattering screen. Both of these values are 
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in stark contrast to the value of 10 ± 1 km s^ (relative to the Solar system barycentre) obtained 
from pulsar timing (Table 1). We note that the scintillation-based velocity depends on a number 
of assumptions about the properties of the effective scattering screen. In contrast, the proper 
motion measurement has a clear and unambiguous timing signature, although the transverse 
velocity itself scales with the pulsar distance. Even allowing that unmodelled effects of Earth 
motion could affect the published scintillation velocities by about 30 km s^, the dispersion- 
based distance would need to be underestimated by a factor of several to make the velocities 
consistent. We believe this is very unlikely, particularly as the tentative detection of a parallax 
gives us some confidence in the dispersion-based distance estimate. Hence, we believe that our 
timing results for both inclination angle and transverse velocity are less susceptible to systematic 
errors and are therefore more secure than those based on scintillation. 

We note that, with the inclination angle being significantly different from 90°, gravitational 
lensing effects (l29t can be neglected. The implied low space velocity, the comparatively low 
derived mass for B and the low orbit eccentricity are all consistent with the idea that the B pulsar 
may have formed by a mechanism different to the usually assumed core-collapse of a helium 
star (OOlU/l . A discussion of its progenitor is presented elsewhere (l32b . We also note that, 
as expected for a double-neutron-star system, there is no evidence for variation in dispersion 
measure as a function of orbital phase. 

Future tests. In contrast to all previous tests of GR, we are now reaching the point with PSR 
J0737— 3037A where expressions of PK parameters to only IPN order may not be sufficient 
anymore for a comparison of theoretical predictions with observations. In particular, we have 
measured uj so precisely (i.e., to a relative precision approaching 10^^) that we expect correc- 
tions at the 2PN level (O to be observationally significant within a few years. These corrections 
include contributions expected from spin-orbit coupling ^33\\34\i . A future determination of the 
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system geometry and the measurement of two other PK parameters at a level of precision sim- 
ilar to that for oj, would allow us to measure the moment of inertia of a neutron star for the 
first time U2\\35l . While this measurement is potentially very difficult, a determination of A's 
moment of inertia to a precision of only 30% would allow us to distinguish between a large 
number of proposed equations of state for dense matter (I36II57I) . The double pulsar would then 
not only provide the best tests of theories of gravity in the strong-field regime as presented here 
but would also give insight into the nature of super-dense matter. 
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with an inset showing an expanded view of the region of principal interest. 
Fig. 1. The tests of general relativity parameter summarized in a graphical form. Constraints on 
the masses of the two stars (A and B) in the PSR J0737— 3039A/B binary system. Shaded re- 
gions are forbidden by the individual mass functions of A and B since sin i must be < 1. Other 
constraining parameters are shown as pairs of lines, where the separation of the lines indicates 
the measurement uncertainty. For the diagonal pair of lines labelled as R, representing the mass 
ratio derived from the measured semi-major axes of the A and B orbits, the measurement pre- 
cision is so good that the line separation only becomes apparent in the enlarged inset, showing 
an expanded view of the region of principal interest. The other constraints shown are based 
on the measured post-Keplerian (PK) parameters interpreted within the framework of general 
relativity. The PK parameter u describes the relativistic precession of the orbit, 7 combines 
gravitational redshift and time dilation, while Pb represents the measured decrease in orbital 
period due to the emission of gravitational waves. The two PK parameters s and r reflect the 
observed Shapiro delay, describing a delay that is added to the pulse arrival times when prop- 
agating through the curved space-time near the companion. The intersection of all line pairs is 
consistent with a single point that corresponds to the masses of A and B. The current uncertain- 
ties in the observed parameters determine the size of this intersection area which is marked in 
blue and which reflects the achieved precision of this test of GR and the mass determination for 
A and B. 

Fig. 2. Measurement of a Shapiro delay demonstrating the curvature of space-time. Timing 
residuals (differences between observed and predicted pulse arrival times) are plotted as a func- 
tion of orbital longitude and illustrate the Shapiro delay for PSR J0737— 3039A. (a) Observed 
timing residuals after a fit of all model parameters given in Table 1 except the Shapiro-delay 
terms r and s which were set to zero and not included in the fit. While a portion of the delay 
is absorbed in an adjustment of the Keplerian parameters, a strong peak at 90° orbital longitude 



18 



remains clearly visible. This is the orbital phase of A's superior conjunction, i.e. when it is 
positioned behind B as viewed from Earth, so that its pulses experience a delay when moving 
through the curved space-time near B. The clear detection of structure in the residuals over the 
whole orbit confirms the detection of the Shapiro delay, which is isolated in (b) by holding all 
parameters to their best-fit values given in Table 1, except the Shapiro delay terms which were 
set to zero. The line shows the predicted delay at the centre of the data span. In both cases, 
residuals were averaged in 1° bins of longitude. 

Fig. 3. Contour plots of the distribution in the plane of the Shapiro-delay parameter Zs = 
— ln(l — s) and the mass of the B pulsar, me. The contours correspond to 68%, 95% and 99% 
confidence limits. 
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Table 1: Parameters for PSR J0737-3039A (A) and PSR J0737-3039B (B). The values were 
derived from pulse timing observations using the DD dTOt and DDS (17^ models of the timing 
analysis program tempo and the Jet Propulsion Laboratory DE405 planetary ephemeris SJE^ . 
Estimated uncertainties, given in parentheses after the values, refer to the least significant digit 
of the tabulated value and are twice the formal l-a values given by tempo. The positional 
parameters are in the DE405 reference frame which is close to that of the International Celestial 
Reference System. Pulsar spin frequencies v = 1/ P are in barycentric dynamical time (TDB) 
units at the timing epoch quoted in Modified Julian Days. The five Keplerian binary parameters 
(Pb,e,u;,To, and x) are derived for pulsar A. The first four of these (with an offset of 180° 
added to uj) and the position parameters were assumed when fitting for B's parameters. Five 
post-Keplerian parameters have now been measured. An independent fit of uj for B yielded a 
value (shown in square brackets) that is consistent with the much more precise result for A. The 
value derived for A was adopted in the final analysis (see (ITib ). The dispersion-based distance 
is based on a model for the interstellar electron density (E?b . 
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Timing parameter 


PSRJ0737-3039A 


PSRJ0737-3039B 


Right Ascension a 


07*^37^51". 24927(3) 




Declination <5 


-30°39'40". 7195(5) 




Proper motion in the RA direction (mas yr~^) 


-3.3(4) 




Proper motion in DecUnation (mas yr"^) 


2.6(5) 




Parallax, tt (mas) 


3(2) 




Spin frequency v (Hz) 


44.054069392744(2) 


0.36056035506(1) 


Spin frequency derivative z> (s~^) 


-3.4156(1) X 10-16 - 


-0.116(1) X 10-15 


Timing Epoch (MJD) 


53156.0 


53156.0 


Dispersion measure DM (cm~*pc) 


48.920(5) 




Orbital period Pj, (day) 


0.10225156248(5) 




Eccentricity e 


0.0877775(9) 




Projected semi-major axis x = (a/c) sini (s) 


1.415032(1) 


1.5161(16) 


Longitude of periastron uj (deg) 


87.0331(8) 


87.0331 + 180.0 


Epoch of periastron To (MJD) 


53155.9074280(2) 




Advance of periastron w (deg/yr) 


16.89947(68) 


[16.96(5)1 


Gravitational redshift parameter 7 (ms) 


0.3856(26) 




Shapiro delay parameter s 


0.99974(-39, +16) 




Shapiro delay parameter r (ps) 


6.21(33) 




Orbital period derivative Pf, 


-1.252(17) X 10-12 




Timing data span (MJD) 


52760 - 53736 


52760 - 53736 


Number of time offsets fitted 


10 


12 


RMS timing residual a (/isec) 


54 


2169 


Total proper motion (mas yr~i) 


4.2(4) 




Distance d(DM) (pc) 


~ 500 




Distance d{n) (pc) 


200 - 1000 


Transverse velocity (d = 500 pc) (km s~^) 


10(1) 




Orbital inclination angle (deg) 


88.69(-76,+50) 


Mass function (Mq) 


0.29096571(87) 


0.3579(11) 



Mass ratio, R 1.0714(11) 
Total system mass (Mq) 2.58708(16) 
Neutron star mass (me) 1.3381(7) 1.2489(7) 
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Table 2: Four independent tests of GR provided by the double pulsar. The second column lists 
the observed PK parameters obtained by fitting a DDS timing model to the data. The third 
column lists the values expected from general relativity given the masses determined from the 
intersection point of the mass ratio R and the periastron advance u;. The last column gives the 
ratio of the observed to expected value for each test. Uncertainties refer to the last quoted digit 
and were determined using Monte Carlo methods. 



PK parameter 



Observed GR expectation Ratio 



7 (ms) 



1.252(17) 1.24787(13) 1.003(14) 

0.3856(26) 0.38418(22) 1.0036(68) 



s 



0.99974(-.39.+ 16) 0.W987(-1S.+ 13) 0.99987(50) 
6.21(33j 6.I53(26J 1.009(55) 
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Fig. 2 
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Supporting Online Material 

1 Observing systems 

The experimental data presented in the main paper are based on pulsar timing observations 
at several frequencies between 320 MHz and 3100 MHz using the Parkes radio telescope in 
Australia, the Lovell radio telescope at Jodrell Bank Observatory, UK, and the Green Bank 
Telescope (GBT) in the USA, between 2003 April and 2006 January. Details of the observing 
systems are summarized in Supporting Table 1 . 

At the Parkes 64-m radio telescope observations were carried out using the centre beam of 
the 20-cm multibeam receiver and a coaxial 10cm/50cm receiver. For each of these cryogeni- 
cally cooled receivers, two orthogonally polarized signals were amplified and down-converted 
to an intermediate frequency. These signals were transferred to band splitters and fed into a 
filterbank system (FB) for each polarization of each feed. The output of each filter was detected 
and summed with its corresponding polarization pair. These summed outputs were high-pass 
filtered and integrated for the sampling interval of 80 fis and then one-bit digitised. While the 
original frequency channels were folded with a reference frequency corresponding to the band 
centre, timing was performed on sub-bands. 

Observations at the GBT utilized two different data acquisition systems. The Berkeley- 
Caltech Pulsar Machine (BCPM) is a flexible filterbank system Q, with which we collected 4- 
bit summed-polarization data. The Green Bank Astronomical Signal Processor (GASP) carries 
out 8-bit Nyquist-sampling of the incoming dual-polarization signal, after which it performs 
coherent dedispersion in software on a Linux-based cluster for each of several 4-MHz channels 
(ElQl). The data stream is then detected, and the two polarizations are usually flux-calibrated 
before summation using a diode noise source as a reference. 

At Jodrell Bank we used a incoherently dedispersing filterbank system. Its parameters are 
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summarized in Table 1, while details of the observing system can be found in ref. 4. 

2 Dedispersion 

Since the interstellar medium (ISM) is ionized, the propagation speed of radio pulses depends 
on their radio frequency with pulses emitted at a high radio frequencies arriving earlier than 
low-frequency pulses. Unless this effect is accounted for, pulses will be broadened over the 
finite observing bandwidth. Two dedispersion techniques are in use. For "incoherent dedisper- 
sion", the bandwidth is sub-divided into a number of frequency channels which are detected 
and sampled independently. Dispersion smearing is thereby reduced to the smearing across an 
individual filterbank channel. The "coherent dedispersion" technique involves the application 
of an inverse "ISM-filter" to the raw voltage data received from the antenna fS). This technique 
is computationally more intensive but removes the effects of dispersion completely. 

At Parkes and Jodrell Bank we obtained incoherently dedispersed data using the filterbank 
systems listed in Table 1. The resulting profiles were summed across frequency channels with 
appropriate delays to remove the effects of interstellar dispersion. For the wide-bandwidth 
Parkes data, where the original frequency channels were folded with a reference frequency cor- 
responding to the band centre, timing was performed on sub-bands. The number of sub-bands 
was chosen such that the dispersion delay across the sub-bands was significantly smaller than 
the overall timing precision. Analysis of TOA data separately for the different sub-bands prop- 
erly accounts for the fact that data at different frequencies received at a given time correspond 
to different orbital phases at emission due to the differential dispersion delay (see e.g. @). 

At the GBT, the BCPM data were divided in four frequency sub-bands, separately dedis- 
persed, folded and timed. In contrast, each GASP 4-MHz channel was coherently dedispersed 
and folded using the channel centre frequency as a reference. The GASP channels were then 
summed appropriately to give a single TOA for each integration. 
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3 Pulse Time-of-Arrival analysis 

Pulse times-of-arrival (TOAs) were computed by correlating the observed pulse profiles with 
synthetic noise-free templates (see Figs. 1 and 2; cf. ref. (7)). All datasets obtained at different 
epochs and frequencies with different data acquisition hardware and telescopes were studied 
for possible systematic errors and artificial correlations. Firstly, correlations between succes- 
sive TOAs were investigated by computing the post-fit root-mean- square (rms) timing residuals 
with averaging of consecutive TOAs, expecting that the rms residual should decrease with the 
square-root of the number of averaged TOAs. Datasets with significant deviations from this 
expected scaling were excluded from the analysis. Secondly, for the GBT observations where 
we recorded data with two different data acquisition systems in parallel, we preferred to use to 
more accurate GASP data and only used BCPM data if no GASP TOAs were available within 2 
minutes of a BCPM TOA. Thirdly, the uncertainties of the TOAs in the remaining datasets were 
studied by inspecting the reduced achieved in the fit of the timing model. For most datasets 
we applied a small quadrature addition and a scaling factor to the uncertainties to obtain the ex- 
pected value of Xred = 1- No adjustments to the TOA uncertainties were needed for the GASP 
data; this is not surprising as the 8-bit sampling provides excellent profile fidelity. Finally, all 
retained datasets were combined in a weighted least-squares fit of the DD and DDS models. 
Following these fits, we verified that the Xred each data subset was still close to unity. A 
total of 131,416 arrival times were included in the final analysis of A while 507 TOAs were 
used for B, most at frequencies close to 820 MHz and 1400 MHz. The much smaller number 
of TOAs for B results from several factors: JBO data were not used, the integration time for B 
was a factor of ten larger than for A, the data were summed over the entire observed frequency 
band, only about 20% of the orbit was used and finally, even in the analysed regions, B was 
often too weak to give a significant TOA. Figures 3 and 4 summarise the TOA distributions for 
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the different observatories for pulsars A and B respectively. Finally, we present the covariance 
matrix as computed by TEMPO for the fit of the DDS timing model in Table 2. 
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Supporting Table 1. Summary of the observing systems used for timing observations of the 
double pulsar. 



Telescope 


Instrument 


Centre 


Gain 


T'sys 


Sample 


Bandwidth 


Number 






freq. (MHz) 


(K/Jy) 


(K) 


interval (/is) 


(MHz) 


of chaimels 






680 


0.66 


45 


80 


64 


128 


Parkes 


FB 


1374 


0.74 


22 


80 


256 


512 






3030 


0.62 


28 


80 


768 


256 


GBT 


BCPM 


820 


2.0 


25 


72 


48 


96 






1400 


2.0 


20 


72 


96 


96 




GASP 


340 


2.0 


70 


0.25 


16 


4 






820 


2.0 


25 


0.25 


64 


16" 






1400 


2.0 


20 


0.25 


64 


16" 


Jodrell Bank 


FB 


610 


1.1 


32 


44.4 


8 


32 




FB 


1396 


1.1 


32 


44.4 


64 


64 



"The number of channels and hence bandwidth that was used varied occasionally within a given session due to the removal of channels 
contaminated with radio frequency interference and/or occasional recording disk-space limitations. 
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Supporting Table 2. Covariance matrix as computed by TEMPO for a fit of the DDS timing 
model to the TOAs of A. 





V 


i> 


Dec 


RA 


PMDec 


PMRA 


X 


e 


To 




w 




7 


DM 


7r 




u 


1.00 



































-0.76 


1.00 






























Dec 


0.16 


-0.31 


1.00 




























RA 


0.10 


-0.08 


0.18 


1.00 


























PMDec 


-0.25 


0.39 


-0.83 


-0.16 


1.00 
























PMRA 


0.01 


-0.28 


0.04 


-0.71 


0.12 


1.00 






















X 


-0.02 


0.02 


0.01 


0.00 


-0.01 


-0.01 


1.00 




















e 


0.00 


0.00 


0.01 


-0.01 


-0.02 


-0.01 


0.66 


1.00 


















2^0 


-0.54 


0.43 


-0.01 


0.02 


0.02 


-0.03 


0.00 


0.01 


1.00 
















Pb 


047 


-0.47 


0.01 


-0.02 


-0.02 


0.02 


-0.15 


-0.15 


-0.85 


1.00 














U) 


-0.54 


0.42 


-0.01 


0.02 


0.02 


-0.03 


0.06 


0.01 


0.99 


-0.84 


1.00 












(jj 


047 


-0.48 


0.01 


-0.02 


-0.02 


0.02 


-0.15 


-0.15 


-0.85 


1.00 


-0.84 


1.00 










1 


-0.02 


-0.02 


0.00 


0.01 


0.00 


0.00 


0.44 


0.01 


-0.03 


0.02 


0.10 


0.03 


1.00 








DM 


-0.02 


0.02 


-0.01 


0.01 


0.03 


0.02 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


1.00 








0.05 


-0.02 


-0.21 


0.16 


0.20 


-0.12 


-0.03 


-0.04 


-0.01 


-0.01 


0.00 


-0.01 


0.02 


-0.02 


1.00 




■Pb 


OOl 


-0.02 


0.00 


0.03 


0.00 


-0.02 


0.05 


0.00 


-0.03 


-0.06 


0.01 


0.01 


0.13 


0.00 


-0.01 


1.00 




-0.01 


0.02 


0.01 


-0.01 


-0.01 


0.00 


0.55 


0.44 


0.01 


-0.10 


0.01 


-0.10 


-0.01 


0.00 


-0.01 


0.00 


1712 


OOl 


-0.03 


-0.01 


0.00 


0.01 


0.01 


-0.85 


-0.75 


-0.02 


0.18 


-0.02 


0.18 


-0.02 


0.00 


0.03 


0.00 



1.00 
-0.72 
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Supporting Figure 1. Pulse profile templates used for TOA determinations for pulsar A. 

Supporting Figure 2. Regions of orbital phase (hatched) used for timing of pulsar B and pulse 

profile templates for these phases derived from and used for the 820 MHz GBT observations 
in May 2005. Each of the template plots covers a range of 60/360 — 0.17 in pulse phase. 
Similar but different templates were used for other frequencies and epochs. While B was clearly 
detectable in these three regions, it is actually brightest in the two cross-hatched regions, but 
because the shape of the profile evolves quickly and dramatically in these regions, they were 
excluded from the timing analysis. 

Supporting Figure 3. Timing residuals obtained for pulsar A for the three telescopes and their 
distribution. The upper panel shows the distribution of observations in frequency. 

Supporting Figure 4. Timing residuals obtained for pulsar B for Parkes and the GBT and their 
distribution. The upper panel shows the distribution of observations in frequency. 
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